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The skeleton formalism aims at extracting and quantifying the filamentary struc- 
ture of the universe is generalized to 3D density fields; a numerical method for com- 
putating a local approximation of the skeleton is presented and validated here on 
Gaussian random fields. It involves solving equation ( H Vp x Vp ) = where Vp and 
Ti. are the gradient and Hessian matrix of the field. This method manages to trace 
well the filamentary structure in 3D fields such as given by numerical simulations of 
the dark matter distribution on large scales and is insensitive to monotonic biasing. 
Two of its characteristics, namely its length and differential length, arc analyzed for 
Gaussian random fields. Its differential length per unit normalized density contrast 
scales like the PDF of the underlying density contrast times the total length times a 
quadratic Edgeworth correction involving the square of the spectral parameter. The 
total length scales like the inverse square smoothing length, with a scaling factor given 
by 0.21(5.28 + n) where n is the power index of the underlying field. This dependency 
implies that the total length can be used to constrain the shape of the underlying 
power spectrum, hence the cosmology. 

Possible applications of the skeleton to galaxy formation and cosmology are discussed. 
As an illustration, the orientation of the spin of dark halos and the orientation of 
the flow near the skeleton is computed for dark matter simulations. The flow is lam- 
inar along the filaments, while spins of dark halos within 500 kpc of the skeleton are 
preferentially orthogonal to the direction of the flow at a level of 25%. 



1 INTRODUCTION 



Recent galaxy surveys like 2dF (CoUess & al.2003) or SDSS 
(Gott & al. 2005) emphasized the complexity of the mat- 
ter distribution in the universe which presents large scale 
structures such as filaments, clusters or walls on the bound- 
aries of low density bubbles (voids) . On the theoretical side, 
the currently favoured scenario suggests that the universe 
evolved from Gaussian initial conditions to form the struc- 
tures that are observed nowadays. Numerical simulations 
have successfully -both statistically and visually, captured 
the main features of the observed filamentary distribution. 

Novikov, Colombi & Dore 2005 (NCD) introduced the 
skeleton formalism in 2D, which aims at extracting and 
analysing the filamentary structure of a given density field. 
This paper extends it to three dimensions in order to de- 
scribe the universe's large scale matter distribution and its 
dynamical environment. 

In the literature, various steps towards a quantitative 
description of the large structures have been suggested. Sta- 
tistical tools such as correlation functions (e.g., Peebles 



1980) and power spectra (e.g. Peacock 1998) have been 
widely used and have been successful in describing mat- 
ter distribution and constraining cosmological parameter. 
Recently, fast algorithms have been designed for first and 
second order (Szapudi & al. 2005), as well as higher order 
statistics (counts in cells etc..) as in (Croton & al. 2004) or 
(Kulkarni & al. 2007). Minkowski functionals (see, e.g., Ker- 
scher 2000 for a review) are so-called "shape finders" (Sahni 
& al. 1998) attempt instead to describe the topology of a 
distribution (see (Sheth & Sahni 2005) for application to 
large scale structure) and can discriminate between Gaus- 
sian and non Gaussian fields as shown in (Doroshkevich & 
al. 1970), (Gott & al. 1986) or more recently (Hikage & al. 
2006)). These topological and statistical estimators analyse 
the distribution of observed galaxies globally and uniformly, 
and make little attempt at recovering the precise geometry 
of the matter distribution, i.e. they do not focus on specific 
regions (such as clumps, voids and filaments). 

Focusing on the identifiable regions of the universe, the 
peak patches theory (Bond & Myers 1996) attempts to de- 
scribe cosmic structures formation through the identifica- 
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tion of the collapse of the dense regions near the density 
peak and surrounding patches. In this framework, the evo- 
lution of patches hierarchy can be understood from the mea- 
surement of only a few characteristics of the patches, while 
assuming that their flow does not depend on their internal 
non-linear dynamics. This line of thought has been extended 
in the Cosmic Web paradigm (Bond, Kofman, Pogosyan, 
1996), which has emphasized that the large scale spatial dis- 
tribution of galaxy clusters and the filaments between them 
can be understood as mildly non-linear enhancement of high 
density peaks and filamentary ridges present in the initial 
gaussian density field. Recently, Hanami (2001) presented 
the so-called skeleton tree formalism: it analyses the pro- 
cess of hierarchical merging and extends the language of the 
peak patch through the analysis of the ridges of the density 
field in an abstract space corresponding to the usual three 
dimensions augmented by the smoothing length. 

The structure of voids in the large scale dark mat- 
ter distribution also has an extended history of theoretical 
modeling (see e.g. Hoffman & Shaham 1982, Icke 1984 or 
Bertschinger 1985) while various void identifiers have been 
designed (see e.g. Platen 2007 and references therein). 

One of the first attempts to develop an algorithm to de- 
tect and trace the filaments in the particle distribution has 
been the minimal spanning tree (MST) technique proposed 
by Doroshkevich in (Doroshkevich & al 1970, Barrow & al 
1985). Starting from a point distribution (a galaxy survey 
or a dark matter simulation), this method constructs the 
graph that connects all the dots with the property of never 
forming closed paths and being of minimal total lengths. 
Interesting statistical features can be extracted from it like 
the shape of the clusters or the length of the trunk (the 
longest path) and branches which are characteristic of the 
filamentarity of the distribution. 

The three-dimensional skeleton described in this paper 
focuses on the critical lines of a distribution, i.e. the set of 
lines joining the critical points in order to be able to com- 
pute the characteristic features of the underlying field (such 
as the total length of the filaments in a cosmological dark 
matter distribution). The skeleton provides a simple math- 
ematical definition of the filaments of a density field based 
on Morse theory (see, e.g., Milnor 1963; Colombi, Pogosyan 
& Souradeep 2000; Jost 2002, Novikov et al. 2006) and thus 
allows their extraction as well as their characterisation. 

Section 2 defines the local skeleton of large scale struc- 
tures. Section 3 introduces the numerical algorithm for con- 
structing the local skeleton, and discusses its properties near 
the critical points (Appendix A gives a more detailed de- 
scription of the algorithm). Section 4 investigates the evo- 
lution of its differential and total length as a function of 
the properties of the underlying field. Appendix C sketches 
the derivation of this differential length. Possible applica- 
tions to cosmology and galaxy formation are discussed in 
section 5, where two illustrations regarding the nature of 
the dark matter fiow near the skeleton are given. 



2 THE LOCAL SKELETON: THEORY 

A comprehensive definition of the skeleton and how its local 
approximation in two dimensions is derived can be found in 



Novikov, Colombi & Dore 2006. To sum up, the so-called 
"real" skeleton is by definition the subset of critical lines 
joining the saddle points of a field to its maxima while fol- 
lowing the gradient's direction (while critical lines link all 
kinds of critical points together). It is easy to picture that 
applying this definition to a 2D field (an altitude map in a 
mountainous region for instance) allows the extraction of the 
ridges of that distribution. Although simple in appearance, 
this definition presents the drawback that it is in essence 
non-local: the presence of the skeleton in a given sub-region 
may depend on the presence of a saddle point in a different 
sub-region. In order to enforce locality, an approximation 
can in fact be derived using Taylor expansion in the vicinity 
of the critical points (i.e. local maxima and saddle points), 
leading to a second order approximation of the skeleton: the 
local skeleton. 



2.1 The 2D local skeleton 

Defining the local critical lines as the set of points where 
the gradient of the field is an extremum along an isodensity 
contour, it can be shown (Novikov, Colombi & Dore 2006) 
that this set of points obeys the equation: 
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where ri and r2 denote space coordinates and p{ri,r2) is 
the density field. Equation ([TJ can be rewritten 



5 = det(HVp,Vp) =0, 



(2) 



where Ti. = d^p/dridr2 is the Hessian (second derivatives 
matrix) of the field. This can be interpreted mathematically 
as the set of points where the gradient of the field is an 
eigenvector of the Hessian (that is, gradient and main 
curvature axis are aligned), which is clearly a local property 
of the field. 

However, in order to correspond to the "real" skeleton 
that traces the ridges and clumps of the field (its structure), 
another condition has to be enforced. For its local approxi- 
mation, it is equivalent to stating that the gradient should 
be minimal (every point of the local skeleton of coordinates r 
should also be a local minimum of the isodensity contour at 
density p (r)). That is, one has to enforce the condition that 
the second eigenvalue of the Hessian should be negative: 

A2<0, and nVp = \iVp, (3) 

where Ai are the eigenvalues of the Hessian and A2 < Ai. 

2.2 The 3D local skeleton 

Let us now derive the generalization of the notion of the 
local skeleton to a three-dimensional space. The philosophy 
is essentially the same but minor differences arise which are 
addressed here. 

Starting from the same definition as in 2D, the skeleton 
should be the set of points where the density is an extremum 
along an isodensity contour. Let {u, v) be a coordinate sys- 
tem along an isocontour (ri {u, v) , r2 {u, v) , ra (u, v)) where 
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Figure 1. Definition of the coordinate system on an isocontour. 

ri,i G {1--3} are the three space coordinates. The definition 
of an isocontour impUes that : 

dp dri dp dr2 , dp drg 
dri du dr2 du 
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= 0. 



dp dri ^ dp dr2 
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Moreover, as the gradient of the field p has to be an 
extremum: 

d 



and 



dv 



(ivpr) = o 



(5) 



Using equations @ and let us derive the equation 
of the local critical lines, which should only depend on 
the field and its first and second order spatial derivatives, 
similarly to equation ([T}. To do so, a coordinate system 
along the isocontour is needed but, as opposed to the 2D 
case, any coordinate system defined on an isocontour will 
be singular in some place as the isocontour is a closed 
surface. In order to avoid this problem, we choose to define 
three coordinates systems and swap from one to another 
when it becomes singular. 

Defining Si three one- dimensional coordinates systems 
so that for different values of Si, one remains in the plane 
(rj, Tk) where i ^ j ^ k and k £ {1..3}. The coordinates 
system Si is singular wherever Vpari. The constrain is to 
satisfy equations. Q and (O for u = Si and v = Sj with 
i ^ j- For any Si, these read: 

dp drs 



dp dn dp dr2 ^.^ , , 

0, and -^-^ + -^-7— + -(6) 
dri dsi dr2 dsi drs dsi 



Choosing i 7^ j 7^ k £ {1..3}, this system becomes after 
some algebra: 
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Indeed, equation ((Tjl reduces equation ([T]) in the 2D case, 
when assuming that the field is constant in the direction 



orthogonal to that 2D plane (the first two terms of equation 
(O are the same as in equation (HJ). The local critical lines 
are thus the set of points that satisfies: 



S = 



= 0, z / j e {1,2,3}. 



(8) 



It is interesting to note that, as in the 2D case, equation 
((8| defines the local critical line as the set of points where 
the gradient of the density is an eigenvector of its Hessian 
matrix (the gradient and the principal curvature axis are 
coUinear) : 



[H-VpxVp) = 0. 



(9) 



Once again, in order to require that the skeleton traces only 
the ridges of the distribution (i.e. the filaments in 3D), re- 
trieving the subset of local critical lines that define the local 
skeleton can be achieved by enforcing a negativity condition 
on the weakest eigenvalues of the Hessian: 



A2<0, A3 < 0, HVp^XiVp 



(10) 



That is, the local skeleton is the subset of the local critical 
where the norm of the 3D gradient is minimal along the 2D 
isodensity contours (as opposed to simply extremal). Note 
that from equation ((SJ it is straightforward to show that any 
monotonic function of the field will have exactly the same 
skeleton as the field itself. 



3 IMPLEMENTATION AND FEATURES 
3.1 Implementation 

Equation (|8]) is at the basis of the numerical implemen- 
tation of the local skeleton determination developed here. 
The details of the algorithm are described in Appendix A, 
while the optimal choice of resolution and smoothing is pre- 
sented in Appendix B. All the computations were performed 
using a specially developed C package: SkelEiQ (Skeleton 
Extractor). This package also includes a fiexible OpenGL 
visualization tool that was used for making the figures in 
this paper. 

Figure [2] presents the skeleton obtained for a density field 
sampled from a numerical simulation of dark matter dis- 
tribution on a 50h,~^Mpc box with 512^ particles using 
GADGET-2 (Springel 2005). The lighter colors represent 
denser regions and the blue skeleton appears to match quite 
well what one could identify as the filaments by eye. Note 
that the skeleton is both a tracer of the topology (it links 
a sub-set of the critical points) and the geometry of its un- 
derlying density. Hence it can be used to compare the ge- 
ometrical and topological properties of various fields, e.g. 
the temperature and the dark matter distribution in hydro- 
dynamical simulations. See also Figure A3 for a graphical 
description of how the local skeleton is drawn. 



3.2 The local skeleton branching properties 

Let us now describe some global branching properties of the 
critical lines and the local skeleton. Important ingredients 



^ Available on request from the authors. 
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Figure 2. The final 3D skeleton derived from a 50 Mpc standard ACDM simulation run with Gadget-2 using 512'^ particles. This result 
is obtained after post treating the skeleton using the method described in Appendix lAl 



of the skeleton are the extrema of the field. Indeed, the 
"real skeleton" is defined as a set of critical lines that 
connect maxima to saddle points. Much of the topological 
behaviour of the skeleton is related to the distribution of 
such extremal points. For the local skeleton described this 
paper, the role of the extrema is similar but the whole set 
of critical lines encompass aditionnal branches linking all 
kind of field extrema together. 

Since the local skeleton is based on a local second 
order approximation of the density field, p, its properties 
can be understood through the properties of the gradient 
Vp and Hessian matrix Ti (p) only. The eigenvalues of Ti 
define the local curvature at any point, thus separating 
space into distinct regions depending on the sign of these 
eigenvalues A^. Within a 3D space, as by definition Aj < \i 
if J > j, there exist four of these regions. Let I be the 
number of negative eigenvalues, then the regions where 
I is equal to 0, 1, 2 and 3. This classification applies to 
critical points of the field in particular, where Vp = 0, 
the maxima (7 = 3) and minima (7 = 0) existing within 
local clumps and voids respectively, while two types of 
saddle points can be distinguished: the filaments type sad- 
dle points (for 7 = 2) and the pancake type ones (for 7 = 1). 

Figure |3] illustrates a second order approximation of 



the density field in the vicinity of the field extrema. The 
total set of critical lines form a fully connected path linking 
all the critical points together and exactly six branches 
pass through each of them in the direction of the three 
eigenvectors of the Hessian. Empirically, it is possible to 
picture the typical behavior of the whole set of critical lines. 
Defining E = {0,1,2,3} and considering a given critical 
point where I = n, ii i < j < k £ E~{n}, this critical point 
C„ is usually linked to three other pairs of critical points 
Ci, Cj and Ck (where I — i, I = j and I — k respectively) 
by critical lines aligned with eigenvectors associated to 
eigenvalues Ai, A2 and A3 respectively at point Cn- Most of 
the time, each of these branches connect to critical point d, 
Cj and Ck along the eigenvector associated with eigenvalue 
Ai, A2 and A3 respectively, evaluated at points d, Cj and 
Ck respectively. In this picture, the critical lines can be seen 
as a fully connected path linking all the different regions 
defined by the sign of the eigenvalues of Ti.. 



The overdense filamentary structure correspond to the 
subset of the critical lines that constitute an approximation 
of the "real" skeleton (i.e. the "ridges" of the distribution). 
This part is the one which links maxima (7 = 3) and 
filamentary saddle points (7 = 2). The typical behaviour 
of such lines is the following: In the immediate vicinity of 
a non-degenerate maximum, two branches of the skeleton 
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(a) > Al > Aa > A3 (/ = 0) (b) Ai > > A2 > A3 (/ = 1) 





(c) Al > A2 > > A3 (/ = 2) (d) Al > A2 > A3 > (/ = 3) 

Figure 3. Illustration of a second-order approximation of the 
density field around a maximum (7 = 0), filament (/ = 1) and 
pancake (/ = 2) saddle point and a minimum (7 = 3). The color 
stands for the density, ranging from purple in low density regions 
to red in high density regions. The axes are the eigenvectors of 
the Hessian, and give the direction of the 6 branches of the local 
critical lines going through these critical points (i.e. where the 
gradient of the field and the eigenvectors of H, are aligned). The 
skeleton is the subset of these critical lines linking maxima (fig. 
|3(a)[ l and filament saddle points (fig. |3(b)| l, in the direction of the 
eigenvector associated with Ai. 



exist, stretching in the eigendirection that corresponds 
to Al. Following one of the branches, denoting as A|| an 
eigenvalue whose eigenvector is parallel to the skeleton 
and Ax, 1,2 as two eigenvalues associated to eigenvectors 
in the perpendicular directions. Near the maximum, 
> Aj| = Al > Ax,i > Ax, 2. As one follows a branch one 
probable outcome is the change of sign of A||, in which 
case the branch will typically end in a saddle point of a 
filamentary type along its Ai direction. There is always 
another branch that starts from this saddle point on the 
other side, thus this type of branches have a fully connected 
structure. However, another possible outcome is that one 
of the orthogonal eigenvalues changes faster than A|| as 
one moves away from the maximum and becomes positive 
before the saddle point is reached. In this case the branch 
of the local skeleton formally terminates, which however in 
reality often means that the skeleton splits at this point in 
two new branches. 

Such branching of the skeleton is especially frequent 
near the maxima of the field, where it accounts for how mul- 
tiple filamentary sections can end up in a single dark matter 
halo. Studiing how skeleton segments merge is relevant for 
questions such as the multipole structure of matter inflow 
onto dark halos (Aubert, Pichon & Colombi 2004, Pichon 
& Aubert 2006). This property of skeleton segments to end 
outside of the critical points is specific to the local defini- 
tion of the skeleton, in contrast to the "real" skeleton whose 
segments are always connected on both ends. 



4 THE SKELETON LENGTH FOR 

SCALE-FREE GAUSSIAN RANDOM FIELDS 

Before considering general cosmological density fields, the 
local skeleton of scale free Gaussian random fields p with null 
average value (p) = will be investigated. For convenience, 
it is useful to define some spectral parameters that depend 
on the spectral index n and on the smoothing length. In the 
statistical description of the skeleton of a random density 
field (Appendix C), the following spectral parameters appear 
to play a role: 
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This introduces three linear scales into the skeleton theory 
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where the first two have a well-known meaning of typical 
separation between zero-crossing of the field _Ro and mean 
distance between extrema, i?« (Bardeen & al. 1986) , while 
the third one, R is, by analogy, the typical distance between 
the inflection points. 

Out of three scales two dimensionless ratios may be con- 
structed that are intrinsic parameters of the theory 
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(16) 



where 7 says how frequent encountering a maximum be- 
tween two zero crossings of the field is, while 7 describes, 
on average, how many inflection points are between two ex- 
trema. For Gaussian fields, these parameters can be easily 
calculated from the power spectrum. Both 7 and 7 range 
from zero up to one. For reference, for the power-law spec- 
tra with index n > —3, smoothed at small scales with a 
Gaussian window, 



/n + 3 . In + h 



(17) 



Note that cosmologically relevant density power spectra 
have n > —3 and thus, while 7 can attain low values, 7 
are always close to unitjQ . 

Appendix C introduces a statistical description of the 
skeleton for the Gaussian and non Gaussian random field. 
This section presents the numerical measurements of the 
properties of the skeleton for scale free Gaussian fields. 

The first quantity of interest is the total length of the 
skeleton, Ltot . In the context of cosmology, Ltot can be linked 
to the total length of the filaments linking clusters together 
and in that sense refiects the history of matter accretion 
as well as the initial distribution of matter (which is sup- 
posed to be similar to a Gaussian random field with a scale- 
dependent effective spectral index similar to the ones consid- 
ered here). Figure |4] presents the result of the measurement 



Cosmological density fields, therefore, have of order one inflec- 
tion point per extremum, unlike, for, example, a mountain range, 
where one encounters many infiection points on a way from a 
mountain top to the bottom 
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Figure 4. Total length L of the skeleton per unit box size for dif- 
ferent smoothing lengths cr = 0.020, 0.027, 0.035; measured over 
25 realisations of Gaussian random fields as a function of the 
spectral index n. While L depends linearly on the spectral index 
n, it grows as a power of a. The dotted lines represent the fits 
obtained using the function: L ^ 0.21(n + 5.28)iT~^. 



of the total length Ltot of the skeleton per unit box size as 
a function of the spectral index and for different smoothing 
lengths a (within the range of validity of the algorithm as 
described in Appendix B). These measurements are carried 
over 25 realisations of scale free 256'^ Gaussian random fields 
as a function of the spectral index n. The sensitivity of the 
skeleton to the value of the spectral index is clear on this 
plot and, if Ltot appears to be a linear function of the spec- 
tral index, it is also clear that it grows as a power law of the 
smoothing length. The dotted lines on figure U shows the 
result of such a fit of the data and seems to work very well. 
A very good approximation of Ltot per unit box size is thus 
given by the function: 



= 0.21(n -f 5.28)o-" 



(18) 



As expected, the exponent of a is measured to be ex- 
actly 2. It can be proved with a simple argument that this 
should be the case for scale free Gaussian fields. In fact, for 
such fields, computing the skeleton over a grid of volume 
and smoothed on a scale a is equivalent to computing the 
skeleton on a grid of volume {al)^ while smoothing on a 
scale (aa) and rescaling the result by a factor 1/a. Because 
of the scale invariance, we also have L (a) = a^^L (aa) and 
so L (a) oc a/a^ = a~^. 

Interestingly, the dependence on the spectral index n 
is close to n -|- 5 which argues for filaments being relatively 
straight between extrema, see Appendix C. A visual exam- 
ination of the filaments confirms this picture. 

Now consider the differential length of the skeleton, 
dL/dr]{rj) where ri = p/ao is the normalized density con- 
trast. This quantity represents the expected length of skele- 
ton that can be measured in a given distribution between 
density contrasts rj and rj + dr). Figure [5] shows the normal- 
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n = 0.219 


0.006 


-0.001 


n = -l 0.212 


0.002 


-0.002 


n = -2 0.206 


-0.005 


-0.008 


0.21 ± 0.005 


0.001 ±0.005 


-0.004 ± 0.003 



Table 1. Measured values of the first three non-null terms in the 
Edgeworth expansion, equation I I19II . for three different values of 
the spectral index n = 0, —1, —2. These results are obtained by 
fitting equation 1191 1 on the data presented in figure |5] on which 
the dotted lines represent the fitted function. The measurements 
show very good agreement, whatever the value of n. 



ized function dL/d-q (rj) as a function of the normalized den- 
sity contrast i] from which was subtracted the probability 
distribution function (PDF) of the field (which, within the 
range of sampling and finite volume effects approximations, 
is a Gaussian function). These values were also averaged 
over 25 realisations of Gaussian fields with spectral index 
71 = 0, —1, —2 sampled on 256^ pixel grids and for a smooth- 
ing length a = 0.027. This value was chosen as a compromise 
between finite volume effect and differentiability of the field 
on a grid discussed in Appendix |B] Considering the error 
bars, it is clear that the value of dL/d-q ijf) is directly linked 
to the spectral index n. 

It is shown in Appendix C that dL/drj (ry) can be written 
using an Edgeworth expansion (see also Novikov, Colombi 
& Dore 2005 for the corresponding proof and fit in 2D): 



dL , , Ltot / 2 /„ 
rf;^W = ^exp(-W2 




C2„7'"^i'2n {v/V2) ,(19) 



where Ltot is the total length of the skeleton, Co = 1 and 
H^n are Hermitte polynomials. Figure 4 demonstrates that 
this expansion also works very well in the 3D case. Remark- 
ably, equation (|19|l does not depend on 7 which again argues 
for the picture of a stiff behaviour of the skeleton for cos- 
mological scale invariant density fields (see Appendix C). 
Table [1] presents the values of the first three coefficients C2n 
obtained by fitting the measurements presented in figure 
[S] (the dotted line of figure [S] are the result of these fits). 
Not only does equation (|19|) allows a very good fit of the 
measured data, but it also appears that only the first order 
term is non-null and the differential length of the skeleton of 
a Gaussian random field with spectral parameter 7 is thus 
given by: 



dL 
dr) 



in) 



Ltot 



exp (-r?V2) (1 + O.2I7'' {rf - l)) 



(20) 



The only non-null coefficients in the expansion are thus 
Co = 1 and C2 = 0.21, to be contrasted to C2 = 0.17 in 
the 2D case. Equation (|20|) can be used as a test of non 
gaussianity like any other topological estimator, such as the 
genus, the PDF etc... as discussed in Novikov et al.(2006), 
since departure from the shape of equation pOf) must appear 
when the skeleton's differential length is computed while the 
underlying field is not Gaussiaijf]. 



of course, given the properties of the skeleton, this won't apply 
if the non Gaussianity involves only a (monotonic) bias 
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Figure 5. Difference between the probability distribution func- 
tion (PDF) of the density field and the normalized differential 
length of the skeleton dL /drj as a, function of the density contrast 
r] = p/ctq. Each curve represents the average value and variance 
of the measured value of dL/dr) over 25 different realisations of 
scale-free Gaussian fields, for different values of the spectral in- 
dex n = 0,-1,-2. The dotted curves represent the estimation 
obtained by fitting data using equation I I19I I (sec table [l] for val- 
ues of the parameters). 



For the matter distribution in the universe, the fila- 
ments are overdense regions along which matter flows. In 
that sense, they are less subject to numerical or observa- 
tional noise and contain most of the information about the 
underlying matter distribution. The skeleton length can thus 
be seen as a method for measuring the power spectrum 
which naturally weights information in different regions ac- 
cording to their importance. 



5 ILLUSTRATION: DYNAMICAL 
ENVIRONMENT OF FILAMENTS 

Drawing the skeleton allows us to pin down the nature of the 
flow around the filaments. Indeed one may roughtly define 
three dynamically distinct regions in large-scale structures: 
voids, clusters and filaments. The first two have been in- 
vestigated in some detail. The filaments represent a fairly 
unexplored venue. Beyond the kinematics (velocity distri- 
bution, spin, etc.), the photometric and spectroscopic prop- 
erties of galaxies (colour, age, metallicity etc.), their mor- 
phology (ellipticals versus spirals, Gini number, Asymme- 
try) or the IGM (gas temperature, WHIM detection, frac- 
tion of gas/metals in the filaments etc.), could also be in- 
vestigated as a function of the distance to, and along the 
filaments. 

In this section, two examples simply illustrate how the 
skeleton can be used to explore the environment of filaments 
in cosmological simulations. 



5.1 Dark matter flow near the skeleton 

Figure |6] displays probability distribution functions (PDF) 
of different characteristics of the dark matter flow along the 
skeleton. In order to understand the correlations between 
the filaments and the velocity field, we computed the PDF 
of its angle relative to the skeleton as a function of its 
intensity (top panels), and the PDF of the angle between 
its largest eigenvector and the skeleton as a function of 
the norm of the corresponding eigenvalue (bottom panels). 
These measurements were achieved by first sampling the 
field characteristics on a grid, averaging particles velocities 
V = (v) and dispersion tensor, AV^^- = {{vi — (vi)){vj — {vj))) 
over each cell, and then computing for each segment the 
distance-weighted average of their PDF. Left and right 
panels yield the resulting PDF computed in a 
and 1000/i~^ Mpc dark matter standard ACDM model 
simulation respectively, at redshift z — and using 512'' 
particles. In both cases, the density and velocity fields 
where sampled on a 512^ pixels grid and smoothed over 
Up = 6 pixels (i.e. s = 1.2h~^ Mpc and s = 12h~^ Mpc 
respectively). The skeleton segments being oriented in the 
direction of increasing density, an angle of 9 = means that 
dark matter is fiowing along the filament in the direction of 
higher density regions. 

The fiows appears to be laminar and its amplitude 
increases with scale: this is expected since on larger scales 
the clusters are more massive, the potential difference 
is larger, hence the flow towards them is faster. Most 
dark matter particles have a mean velocity of about 300 
(resp. 400) km/s along the fllament and a dispersion of 
about 100 (resp. 150) km/s orthogonal to the filaments 
for the two scales considered here. The angular spread 
(panels 6(a)|6(by I also increases with scale, from about 
30° to about 45°, reflecting the larger internal heat of the 
filament, also seen in (panels [6(cj|6(d) I . 



The qualitative shape of this PDF may be explained 
by the advection of new halos onto the "highways" cor- 
responding to the mean flow. The first eigenvector of the 
dispersion tensor is on average clearly orthogonal to the 
filament, reflecting the velocity of dark matter falling onto 
the filaments. Note that the distribution is decreasing 
monotonously with in panel 6(a) some dark matter 



particles statistically even move downhill, and their relative 
fraction decreases with scale. The filaments are collecting 
matter away from the underdense regions. Smaller filaments 
empty smaller voids, which tend to get depleted earlier than 
larger ones; hence this may explain why the flow becomes 
more orderly at smaller scale as accretion diminishes. 
Note that the redshift evolution (not shown here) of this 
distribution follows closely its scale evolution, the z — 15 
PDF over lOO/i"^ Mpc resembling the z = PDF over 
1000/i~^ Mpc dSousbie, PhD Thesis, 20061 ). 

The detailed nature of the flow should eventually be 
investigated in a smoothing scale independent manner, in 
order to derive universal features which would only depend 
on the cosmology and the initial power spectrum. Its evo- 
lution with redshift or with the cosmology should also be 
systematically analysed. 
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Figure 6. Top panels: Probability distribution function (PDF) of the velocity field V of the dark matter along the skeleton as a 
function of its angle 8 with the skeleton and its norm. The measurements were achieved on a lOO/i^^ Mpc and lOOOh"^ Mpc dark matter 
simulation featuring 512^ particles and a standard ACDM model, smoothed over a scale s = 1.2h~^ Mpc and s = 12h~^ Mpc (left and 
right panels respectively). The skeleton is oriented in the direction of increasing density. Dark matter appears to be flowing along the 
filaments in the direction of higher density regions (i.e. halos). Bottom panels: PDF of main eigenvector of the velocity dispersion tensor 
AVij as a function of its angle 6 with the skeleton and its eigenvalue amplitude. The peak of the PDF corresponds to high velocity 
dispersion orthogonal to the filaments, which is coherent with the picture of dark matter being accreted orthogonally by the filaments 
before flowing along them. Note the increase in velocity dispersion with scale [left and right panels) as well as the larger angular dispersion 
in the dark matter flow. This trend is also found while considering the same simulation at higher z. 



5.2 Dark matter spin-skeleton connection 

The geometric orientation of the spin of dark matter halos 
corresponds to anotlier feature of tlie large scale structure 
which can be ciiaracterized using the skeleton. The spin of 
dark halos was computed using the classical friend-of-friend 
(FOF) algorithm with 0.2 times the interparticular distance 
as linking length and retaining only halos containing more 
than 100 particles. Figure [7] displays the excess probability 



of alignment of the halos' spins with the closest skeleton 
segment for different distances [0, 0.5], [0.5, 1.5], [1.5, 2.5] and 
[2.5, 3.5]h~^ Mpc. This probability reaches 25 % for an angle 
6 — ■k/2 between the spin and the skeleton: the spin of dark 
matter halos is preferentially orthogonal to the filament they 
belong to. This trend accounts for the fact that the filaments 
are the locus of laminar fiow where halos coalesce along the 
direction of the filaments parallel to the mean flow, hence 
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Figure 7. Excess probability of spin alignment with the 
local skeleton computed from the average of three 512^ 
bOh~^ Mpc ACDM simulations at different distances: d £ 
[0, 500], [500, 1500], [1500, 2500] and [2500, 3500]/i-l kpc from the 
closest skeleton segments. This figure demonstrates that on aver- 
age the spin of dark matter halos tends to be orthogonal to the 
local filaments at a level of 25 % for distances shorter than 500 
kpc. The simulation is analysed at redshift zero. 



acquiring momentum orthogonal to the flow, as observed in 
(Aubert, Pichon & Colombi 2004). 



6 CONCLUSION & PERSPECTIVES 

The 3D skeleton formalism is a weU-defined framework 
for studying the filamentary structure of a distribution. 
The "real" skeleton is defined as the set of critical lines 
joining saddle points to maxima of the field along the 
gradient. A local approximation of it was introduced in 
section [2] along with a numerical method allowing a fast 
retrieval of the locus of the filaments from a sampled field 
(see also appendix This method involves computing the 
null isodensity surfaces of each component of a function 
S — {H ■ Vp X Vp) of the gradient, Vp, and Hessian 
matrix, Ti of this field. 

The ability to localize and characterize the filamentary 
structure of matter distribution in the universe opens 
the prospect of many applications for the skeleton as 
discussed in Section 4 and 5. It has been shown in section 
|4] that for a Gaussian random field, the total length of 
the skeleton per unit volume depended only of the chosen 
smoothing length a and spectral index n, with a specific 
functional from which was both fitted from simulations 
and motivated in Appendix C. In this sense, the local 
skeleton provides a direct measurement of the local shape 
of the power spectrum, P{k), on various scales depending 
on the smoothing applied to the underlying field. Though 



there exist other ways to measure the power spectrum of 
a given distribution, the skeleton length is promising as it 
relies only on the filamentary structure of the distribution. 
The analysis of the length of the skeleton of the galaxy 
distribution in the SDSS as a measurement of cosmological 
parameter fim can be found in HSousbie et al. 2007ap . 

The skeleton may also be used as an isotropy probe. 
It corresponds in fact to a good candidate for the Alcock- 
Paczynski ( |Alcock fc Paczynski 1979| test, since the 
apparent longitudinal to transverse length of skeleton 
segments should directly constrain the curvature of space 
in a manner which is bias-independent. This test will be 
presented in a forthcoming letter (jPichon et al. 2007bp . 

It was demonstrated in section 4 that the dark matter 
flow in the vicinity of filaments was dominantly laminar 
along the filaments and shows signs of orthogonal accretion 
corresponding to the infall of dark matter collected from 
the voids. It also showed that the spin of dark matter halos 
were preferentially orthogonal to the filament's direction, a 
feature which can be understood as a consequence of merger 
events taking place along these filaments. A clear virtue of 
the local skeleton is that since it relies on a local expansion 
of the field, it can deal with truncated/masked fields, 
segmented or vanishing ridges or isolated structures. Note 
finally that the fit, equation I2UI opens the prospect of using 
the local skeleton to estimate the bias in observed surveys. 
The idea is to compute the PDF of galaxies on the one hand, 
which depends on the mass to light ratio of the sample, and 
the differential length (equation pOjl 'l on the other hand. 
Since the former depends on the bias, whereas the later does 
not, comparing the two should give an estimate of the bias. 
On the other hand, the local formulation of the skeleton 
presents some limitations. Mainly, it is not fully connected: 
it has by construction (since it is drawn from a second 
order Taylor expansion of the field) only 2 segments per 
maxima whereas full connection would require 3 or more. 
A consequence is that it cannot represent merging filaments. 

One could also use the curvature and torsion of fila- 
ments as cosmological probes, since the acceleration of the 
universe induced by the cosmological constant is likely to 
straighten the filaments, though the fact that the local skele- 
ton has only two segments near its maxima (the other seg- 
ments must branch out) is likely to introduce some artifacts. 
The topology and geometry of the skeleton near the density 
peaks and the redshift evolution of the skeleton of the large 
scale structures may prove of interest, for instance to study 
the frequency of reconnection, though again the local skele- 
ton is not ideal in this respect. It would also be interesting 
to construct the skeleton in higher dimensions, for instance 
in space-time, to trace the events lines, but again connection 
is critical. In a forthcoming paper, an alternative algorithm 
for the indentification of the skeleton, loosely based on a 
least action formulation, will be presented. It is complemen- 
tary to the solution presented in this paper and will allow 
us to tackle those points for which the local skeleton is less 
efficient. Finally, the 3D skeleton algorithm could possibly 
be applied to other fields of research, such as neurology, in 
order to trace the neural network. 
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APPENDIX A: NUMERICAL 
IMPLEMENTATION 

All the computations were performed using a specially 
developed C package: SkelEjjj (Skeleton Extractor). This 
package also includes a flexible OpenGL visualization tool 
that was used for making the figures in this paper. 

The first step before computing the skeleton requires 
obtaining a density field from a discrete point-like distri- 
bution. This is achieved by smoothing appropriately the 
density field on a grid so that it is not singular (i.e. is 
sufficiently differentiable) but still contains all the topo- 
logical information. The density field is computed using 
Cloud-In-Cell (CIC) interpolation (e.g., R.W. Hockney 
1988) and convolving the result with Gaussian windows of 
different widths. As was shown in section [B] the grid size 
and smoothing length are decisive parameters. It is then 
necessary to compute first and second derivatives of the 
field on the grid, which can be done using finite difference 
or Fourier transform method, giving very similar results. 
In fact choosing one method or the other does not seem to 
have any influence on the resulting skeleton if the field is 
smooth enough (which is anyway a necessary condition). 

The next step involves solving the system of equa- 
tions (|8}; the solution of this system corresponding to the 
intersection of two of the three solutions of equations ([7]). 
This is done by computing the 3D meshes of the two- 
dimensional surfaces that are solution to these equations: 
the skeleton is at the intersection of two of them, depending 
on the value of the gradient at the point considered. Solving 
equation ^ is equivalent to finding the null isocontour 
of field Si, which can be done using the marching cube 
algorithm ( [Lorensen, W. and Harvey E. (1987) 1. The basic 
idea is to consider every cell of the grid as an individual 
cube. One can then compute the value of every Si for the 
eight vertices and it is easy to check whether the isosurface 
intersects the cube or not. In fact, every vertex is above 
or below a threshold value (in this case 0), which gives a 
total of 2* — 256 types of intersections (only 15 of them 
being instrinsically different) that can be precomputed 
as illustrated in figure lAll The exact positions of the 
intersections are computed using quadratic interpolation. 
This yields the position of the intersections of the grid and 
the isocontour, and defines triangles that smartly link those 
intersection vertices: one can then reconstruct a very good 
approximation of what the isocontour is. 

Which surfaces should be used for each cell is decided 
by computing dk = det (ri, rj, Vp) , i ^ j k £ {1,2,3} 
and selecting only the two Sk for which dk is maximal. This 
gives two surfaces defined by triangles whose intersection 
can be efficiently computed: it amounts to computing the 
intersection of triangle pairs only. It is then straightforward 



* Available on request from the authors. 
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Figure Al. Illustration of the different possible configura- 
tions of a grid cell used for marching cube algorithm. Given 
a field / and isocontour / = 0, a blue ball represents a 
vertex where / > 0. It is then easy to build the isocon- 
tour by linearly interpolating the value of / along the edges. 
This picture was borrowed from James Sharman's website, 
|http: / /www.exafiop.org/docs/marchcubes/ind.html| 

to compute the eigenvalue of the Hessian for every segment 
and keep or reject them depending on the previously defined 
criteria (equation HOjl) in order to draw the local skeleton. 
The exact same method was used for efficiently and consis- 
tently finding the extrema and saddle points of the field. In 
fact, if one defines three fields /; = dp/dri, those critical 
points are the intersections of the three isocontour surfaces 
fi = 0. One can then decide if a critical point is a maximum, 
minimum or saddle point by checking the value of eigen- 
values of the Hessian (i.e. curvature). Although marching 
cubes algorithm are very efficient for computing isodensity 
contour, they present some drawbacks for ambiguous con- 
figurations. Indeed, as illustrated on figure IA2I some con- 
figurations are degenerate and one cannot decide where the 
isosurface should pass. This problem happens most of the 
time around critical points where the value of the field can 
go above and below the threshold within one cell. It induces 
the loss a small skeleton segments. 

In order to obtain a smooth skeleton that does not 
present holes and to retrieve connectivity information (i.e. 
to be able to follow the skeleton from one point to another), 
a three steps post-processing is applied. There is of course 
some arbitrariness in this process: here the algorithm is 
based on a weighted marking system to achieve this result 
(where the weights are assigned depending on the relative 
importance of the selection criteria): (i) the branches that 
were missed around the extrema are regenerated using 
the fact that the skeleton around an extremum is along 
the main curvature axis (i.e. along the first eigenvector 
of the Hessian). So for each extremum, marks are given 
to all skeleton segment, favoring those at small distances 
and with similar orientation as the main eigenvector of Ti. 
Each extremum is eventually connected to the segment 
with the highest mark, (ii) the gaps between segments in 
the sequence of skeleton branches are filled. Starting from 
segments connected to extrema, all segments are visited 
iteratively: for the running segment, a mark is now assigned 
to all other unprocessed segments, based upon their relative 
distance, their relative angle, and relative orientation. Note 
that the corresponding cost functions are non hnear: for 




Figure A2. Illustration of a drawback of marching cubes algo- 
rithm. The green surface is an isosurface solution of equation (0 
and the light blue line is the resulting skeleton. The red diamond 
represent a field maximum. It is clear on this picture that the 
algorithm misses the part close to that maximum, thus creating 
a spurious hole in the skeleton. 



instance segments with too large a relative angle are given 
an exponentially negative mark, (iii) finally, all segments 
which have not been considered during step (ii) are dropped. 
The resulting skeleton is shown on figure [2] A detailed 
accounting of all stages of the skeleton extraction, including 
the post treatment is given in ( |Sousbie, PhD Thesis, 2006[ ) 
(which gives the exact marking scheme described above), 
while the code is available upon request from the authors. 

From a performance point of view, this method presents 
the advantage of being both fast and robust. The computa- 
tional cost in fact mainly scales as the number of pixels in 
grid Ng\ the cost of computing the iso surfaces intersec- 
tions is neglectible given the possibility of computing only 
the intersections of faces belonging to the same pixel. It is 
moreover memory-efficient and can be trivially parallelized: 
the computation can be done on sub grid regions and then 
merged. On a modern computer, the memory requirement 
corresponds to the requirement to store one sub-grid and 
its three isosurface, which can be arbitrarily small, and the 
computational time for a 128^ pixels grid is of the order of 
a few seconds on a modern desktop computer while only a 
few tens of minutes is necessary for a grid of 1024^ pixels. 



APPENDIX B: SMOOTHING LENGTH AND 
RESOLUTION 

One aspect of the numerical implementation that deserves 
special attention is the issue of smoothing. In the main 
text, we consider the total skeleton of Gaussian random 
fields, focusing mainly on two of its properties: its length 
L and differential length dL/drj. The algorithm presented 
in this paper deals with the numerical computation of the 
skeleton of a discretized realization of a given field. It is 
thus important in the first instance to be able to deal with 
the influence of this discretization on the measured skeleton 
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(a) The three fields iso-surfaces the intersection of which con- 
stitutes the critical lines 




(c) The local critical lines obtained by selecting only the two 
least degenerate fields depending on the value of the gradient. 




(b) The resulting critical lines made of all the intersection of 
any two of the three iso-surfaces shown in |3(a)| 




(d) The skeleton obtained after enforcing condition llOl on the 
local critical lines: Ai > and A2,A3 < 0. 



Figure A3. Illustration of the process of the skeleton computation. White points are dark matter particles extracted from a standard 
ACDM simulation run using Gadget-2. The skeleton is defined as the intersection of two (among three) iso-surfaces (fig. |3(a)| l. Defining 
the curvature as A^ with TiVp = A^Vp and Vj > i, Aj < A^ (p being the density and H its Hessian), it is possible to select only some 
parts of the skeleton depending on the value of A; and retrieve only the filaments (fig. |3(d)[ l. Using a simple post treatment, it is then 
possible to remove insignificant pieces and obtain the precise locus of the filaments (fig. O- 



properties (see e.g. Colombi & al. 2000). 

The statistical properties of a scale-free Gaussian 
random field can be described using only two numbers; 
its spectral index n and the amplitude A of its power 
spectrum P (k) = Ak^ , where k is the wave number. The 
skeleton formalism is totally independent of the amplitude 
of the field, so only the value of n is of interest to us. 
Consider a realization of a 3D scale-free Gaussian random 
field with spectral index n on a A'^ pixel grid. In order to 
ensure sufficient differentiability, this field is convolved to a 
Gaussian kernel whose scale a is expressed per unit box size. 
The value of a limits the size of the smallest scale that can 
be considered, while the finite size of the grid imposes an 
upper limit. Figure iBll presents the measured value of the 



spectral parameter 7^ = {n + 3)/{n + 5) as a function of a, 
for 25 realizations of Gaussian random fields with spectral 
index n £ {0, —1, —2}, together with the theoretical value, 
measured by integrating the power spectrum truncated to 
grid limit frequencies. As expected, a departure from theory 
is observed for higher values of a, especially for fields with 
lower spectral index where most of the power is concen- 
trated on small values of k (i.e. on large scales). This sets an 
upper limit on the value of the smoothing scale and so we 
will only be considering fields smoothed on scales a < 0.035. 

The other constraint on the value of a arises from the 
fact that the skeleton computation algorithm requires a field 
that is continuously differentiable two times in the finite dif- 
ference scheme sense. This means that the smoothing length 
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should be large enough for the computational errors on field 
derivatives to be neglectible. These considerations imply a 
lower limit on the smoothing length value expressed in num- 
ber of pixels CTp = oNg. In order to estimate this limit, we 
again generated Gaussian random fields with different spec- 
tral indices over a 256^ (A'^g = 256) pixels grid and downsam- 
pled them on grids with eight different values of Ng ranging 
from Ng — 64 up to Ng = 224. Figure iBll presents the evo- 
lution of the measured skeleton length for these realisations, 
each of them being computed for a smoothing scale corre- 
sponding to a constant fraction of the box size a ~ 0.031 
but to different values of ap ranging from o-p = 1 up to 
(jp — 7. One would clearly expect the length of the skeleton 
to depend only on the value of a as long as the numeri- 
cal approximations are neglectible, which seems to be the 
case only for values of ap at least of order 5 pixels. For a 
given sampling Ng, this limits the possible smoothing scale 
to a > 5/Ng. As was noted previously, this exact value de- 
pends on the considered spectral index, so we chose to con- 
sider the worst case, n = 1, where the fluctuations of the 
held do not dampen on small scales thus making the field 
naturally not smooth on any scale. 

In this paper, all fields considered are sampled over Ng = 
256 cubic grids, so in order to respect the constraints de- 
scribed above, the fields are smoothed on scales in the range 
0.02 <a < 0.035. 



Figure Bl. Top: evolution of the measured spectral parameter 
7 (see Eq. I16I I for 25 realisations of Gaussian random fields witli 
spectral index n = (red crosses), n = — 1 (green triangles) and 
n = — 1 (blue discs) as a function of the smoothing scale cr ex- 
pressed in box size units. The three continuous lines represent the 
expected theoretical values, measured by integrating the power 
spectrum truncated to grid limit frequencies. The dotted lines 
are the theoretical expectations (equation (1171 ) without account- 
ing for finite volume effects. For higher values of cr, the finite box 
size effects have more influence and the measured value of 7 tends 
to differ from the correct one, thus limiting the maximal smooth- 
ing scale Bottom: evolution of the measured length of the total 
skeleton in box size units as a function of the smoothing length 
in pixels ap , for different values of the spectral index n and while 
keeping the smoothing scale to a constant fraction of the box size 
cr fa 0.031 . The measurements are obtained by resampling one 
initial realisation of a Gaussian random field (genrerated over a 
256^ pixels grid) on smaller resolution grids and smoothing the 
resulting fields over the appropriate number of pixels. The mea- 
sured length of the total skeleton appears to become stable for 



responds to cr > 0.19 for a field sampled on a 256 pixels grid. 
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APPENDIX C: THE THEORETICAL DIFFERENTIAL LENGTH OF THE SKELETON 
CI The length of the skeleton 

To estimate the length, /I(pth), of the local skeleton per unit volum^fl consider the vicinity of the points through which the 
local critical line passes, 5"; = 0, = 0. (Since the sets of conditions {Si,Sj) = (0,0), i 7^ j is degenerate, without loss of 
generality one can assume a particular choice for i and j). Define the set of points, £, in the excursion p > pth near the critical 
line solutions that satisfy — A5i/2 < Si < ASi/2 and —ASj/2 < Sj < ASj/2 where ASi and ASj are sufficiently small so 
that the linear expansion ASi ~ \/Si ■ dr, ASj ~ VS'j • dr holdfl . The fraction of the total volume the set E occupies (the 
filling factor) is 

V(pth,A5,,A5j)= / dp dS^ dS, d^'iVSi) d''{\/Sj)V{p,S„Sj,VS„VSj), (CI) 

■^P>Pth •^-ASi/2 J-ASj/2 J 

where ^(p, 5;, 5j, V5i, V5j) is the joint probability distribution function of the quantities {p,Si,Sj,\7Si,\7Sj). Here the 
seemingly redundant distribution of the gradients VSi and VSj was introduced to have the expression for the fraction of the 
total volume occupied by a differential subset of £ that has specific values of the gradients \7Si,\7Sj (within d^{\7Si) and 
d'(V5,)) 

dV(pth,A5„A5,,V5„V5,) =d='(V50 rf'(V5,) / dp ds J dSjP{p,S„Sj,VS„VSj) . (C2) 

Jp>Pth J-ASi/2 J-ASj/2 

Since the area, E, of a section locally orthogonal to the such subset, is simply (modulo some trigonometry) given by 
E(A5i, A5j, V5i, V5j) = A5,A5j/|V5, x VSj\, 

dividing dV by E, integrating over all possible gradients \/Si,\7Sj and then taking the limit {ASi, ASj) — > (0,0) yields for 
the skeleton length per unit volume: 



'C(pth) = lim 



dV(pth, A5„ A5,, V5„ V5,) 



(ASi,ASj)^(o,o) J E(A5i, A5j, V5i, V5j) 

dp / d^(V5,) d^iVSj) |V5, X VSj\P{p,S, = 0,Sj = 0, V5,, V5j). (C3) 

p>pth 

This generalizes the calculation of NCD to 3 dimensions; the length of the local skeleton is defined by the knowledge of the 
density field and its partial derivatives up to third order, as expected. In order to understand the scaling involved in the 
computation of £., let us rewrite this equation in terms of dimensionless quantities; 

- - - ^'^P _ d^P 2 _ c 2v7 _V7C II^A\ 

ori OTiOrj ariOTjOrk 

with, following ( |Bardeen, et al., 1986[ ), 

k^dk 



27r2 



P(A:)r", (C5) 



where P{k) is the power-spectrum of p. Equation (|9} and its gradient can be written more conveniently, using the totally 
antisymmetric tensor, e'-'*, as 

Si = ^ e'^'^aJjiXiXfe , and VmSi = V Si{Xk,XkUXklm) = '^^^''^'' [XjlrnXlXk + ^[XjlXlrnXk + XjlXkmXl]) . (C6) 



jkl jkl 



Indeed, expressions for Si and Vsi involve up to third derivatives of the field x. All the quantities defined by Eq. HC4|I are 
dimensionless and their variance do not depend on spectral parameters (they are pure numbers) except Vsi. Keeping that in 
mind, with these new notations, equation (IC3|) becomes 



Cipth)=(—)l dx d\Vai) d^{ysj)\Vs^xVsj\V{x,s^ = 0,Sj =G,Vs^,Vaj). (C7) 



Equation (|C7P is the formal expression for the length per unit volume of the total set of critical lines. 

Let us express the joint distribution function 7^(77, Si, Sj, Vsi, Vsj) in terms of the joint distribution function of the 
underlying field and its derivatives P{x,Xk,Xkm, ■ ■ ■)■ Introducing the 20 components vector composed of x and its successive 
partial derivatives up to third order (see equation IC4|) . or in dimensionless units, X = {x, Xk, Xki, Xkim) (symmetries in the 



the distinction is made here between the theoretical expectation, £(pth)i in this section and the estimator, L in the main text. 
^ In such small neighborhood of a critical line there are no other critical lines. Note, that the linear expansion of Si breaks near the 
extrema of the field, where V5i = 0, which allows several critical lines to intersect at such points. However, extremal points are of 
measure zero as far as the computation of the length of the skeleton is concerned. 
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derivative tensors of second and third order are assumed to be exploited to reduce the effective number of variables), we can 
find "P as a marginalization over the distribution of the field values: 

dxd'^Xk(fxud^''xklmP{x,Xk,Xkl,Xklm)SD{x - ri)^ D {Si(Xk , X kl) - Si)Sn{Sj{Xk,Xkl) - Sj) X 

5D{ySt{Xk,Xkl,Xklm) ~ y Si)SD{ySj{Xk,Xkl,Xklm) - VSj) , (C8) 



which yields the appropriate 9D probability distribution. Putting equation HC8|) into equation HC7P and differentiating with 
respect to 77 = p/o"o = x, and accounting for the two Delta functions in Vsi and Vs-, yields, using equation (|f 6p to rewrite 
01, 1 02 in terms of R: 

dC / f \ ^ [ 

= (^-^j / d:''Xkd^Xkld'''^Xklm\ySiXVSj\P{ri,Xk,XkUXklrn)5D{Si{Xk,Xkl))SD{Sj{xk,Xkl))- (C9) 

Equation (|C9|I is the formal expression for the differential length per unit volume of the total set of critical lines, described 
in the main text. Given its dimensionality, it remains a daunting task to compute the 19D integral in equation HC9|) . Note 
that Vsi and Vsj are now functions of {xk,Xki,Xkim) and Si and Sj are function of {xk,Xki) given by equation HC6p . In 
equation (|C9|I . the two Delta functions couple the differents Xk,Xki while accounting for the fact that the integral should be 
restricted to the intersection of the two iso surfaces, i.e. along the critical lines. The modulus in |Vsi x Vsj| reflects the fact 
that the summation of skeleton segments is not algebraic, which complicates also the reduction of equation (|C9|) . For the 
set of local critical lines, there are no restriction to the region of integration. If one is interested in the local skeleton, the 
integration should be restricted to regions where the condition given by equation (|1U|I holds. 
The total length of the critical lines is 

drj— = J d'^Xkd^Xkld^^Xklm \S/Si XS/Sj\P{Xk,Xkl,Xklm)SDiSi{xk,Xkl))SD{Sj{Xk,Xkl)) ■ (CIO) 



C2 £{pth) for Gaussian random field 

Since a Gaussian field does not correlate with its derivatives of odd orders (this is easy to understand using symmetries in 
Fourier space), the joint distribution function P{x,Xk,Xki,Xkim) can be written as 

P{x,Xk,Xkl,Xkhn) = Po{x,Xkl)Pl{xk,Xklm)- (CU) 

In Po, the only dependence on the power spectrum of the underlying field is in the parameter 7 (c.f. equation (|16p ) that 
describes the correlation between the field and its second derivatives. Similarly Pi(xi,Xijk) only involves 7 which describes 
the correlation between the gradient of the field and its third derivatives. Therefore dC/dr) depends only on 77, 7? 7 and 7, 
as argued in the main text. Note that by symmetry, dC/dr) should then be an even function of 77 for the total set of critical 
lines. The total length of the skeleton, Ctot, which follows from marginalization of the equation (|C9|I over 77 may depend only 
on 7 and R since the integration of Po{r],Xki) over 77 cancels the dependancy over 7. 



C2.1 The "stiff" filament approximation 

The 1/R^ scaling in equation (|C9|I reflects the basic fact that by definition the local skeleton is almost straight within a 
volume that has one inflection point ~ 7£^. A straight segment through such volume has length ~ R, thus the expected length 
per unit volume is ~ 1/R? . The dependence on the spectral index is then 1/7?^ oc (n + 7)/a"^, recalling that a is the smoothing 
length in units of the total boxsize. Is this the behaviour with n that one should expect in simulations ? 
Let us write formally 

VSi X VSj = A(x-fc, Xkl,Xklm) + 76(2;*;, Xkl,Xkhn) + J^C{xk,Xkl) ■ (C12) 

Suppose the last term dominates statisticalljQ. Then, since ^/R — 1/R*, and given that C{xk,xi,n) does not depend on the 
third derivative of the field (which can then be integrated out), equation HC12|) becomes 



dri \rJ J 



d^Xkd'^Xki |C(a;fc,a;;,„)| Po{r),Xki) Pi{xk)5D{,Si{xk,Xki))SD{8j{xk,Xki)) ■ (C13) 



It is easy to foresee when this regime is valid. The same argument as before implies that the scaling arises when the 

skeleton is almost straight within a volume that contains one extremum, This is supported by the fact that the integral 
term does not depend on the third derivatives, thus inflection points play no role, and any dependence on 7 drops out. This 
picture corresponds to a skeleton connecting extrema with relatively straight segments. The scaling is then l/i?^ oc (n + 5)/o"^. 
For the total length of the critical lines, integration over 77 gives 



or equivalently assume that the amplitude of derivative of the Hessian is negligible relative to the amplitude of the Hessian 



© 0000 RAS, MNRAS 000, 000-000 



16 T. Sousbie, C. Pichon, S. Colombi, D. Novikov & D. Pogosyan 



Ctot~(^-^^ J d^Xkd'^Xki \C{xk xim)\ Po{xki)Pi{xk)SD{si{xk,Xki))5n{sj{xk,Xki)) oc {n + 5)a ^, (C14) 

since the integral is just a pure number. This is very close to the scaling with n that was found in the numerical fit, 
equation (|f 8[). The differential length in the stiff regime is then only the function of 7 times Ctot- The upshot of this paragraph 
is to demonstrate the theoretical consistency between the scaling in (n + 5)a"~^ of Ctot and the fact that dC/drj does not 
depend on 7 for scale free Gaussian random fields. 



C2.2 Joint distribution of the field and its derivatives for a Gaussian random field 
The full expression Po{x,Xki) is given in ( |Bardeen, et al., 1986[ ). Introducing variables 



-Aa; = —{xii + X22 + X33) 



V = 2 (^33 — a::ii) , w ■ 



(2X22 — Xll — X33) , 



(C15) 



in place of diagonal elements of the Hessian (xii, X22, 133) one finds that u, 11, ixj, X12, Xia, 123 are uncorrelated. Importantly, 
the field, x is only correlated with it = Aa; and 



(xu) = 7, {xv) = 0, {xw) — 0, {xxki) — 0, k ^ I, 

where 7 is the same quantity as in equation (|16|) . The full expression of Po{x,Xki) is then 



Po{x,Xki)dxd Xki = 



(15) 



5/2 



(27r)V2(l_^2)i/; 



■ exp 



{x — 'yu)^ 
2(1-7^) 



cxp 



15 / 2 2 2 2 2 \ 
- — (?; +w + a;i2 + a;i3 + 2:23) 



(C16) 



dxdudvdwdxi2dxi3dx23, 



and is described by only one correlation parameter 7. 

A similar procedure can be performed for the joint probability of the first and third derivatives of the fields, Pi{xi, Xijk) 
by defining the following nine parameters (see also Hanami 99): 



Ui = ViU, Vi = -e'^ Vi (VjVj — Vk'Vk)x , with j < k, and Wi 



(C17) 



and replacing the variables {xiii,Xi22,Xi33) with (ui,Vi,Wi). In that case, the only cross-correlations in the vector 
{x\,X2,X3, u\,v\,w\,U2,V2,W2,U3,V3,W3,x\23) wMch do uot vanish are between the same components of the gradient and the 
gradient of the Laplacian of the field: 



{xiUi) = 7/3, i = 1,2, 3, 

where 7 is the same quantity as in equation (|16p . This allows us to write: 



p / ^w3 ^10 _ 105^''^3^d^w» d'^Vj dxi23 

l^l(Xi,Xijk)a Xid Xijk- (-27r)13/2(]^ _ .::j,2-)3/2 



105 
2 



2:123 + ^^(vf + w'i) 



exp 



(CIS) 



3(tii — ^XiY 'ix} 
2(1 -72) 2" 



C2.3 The differential length for Gaussian fields 



What is the dependence of the skeleton differential length on the parameter 7 and the threshold rjl Let us look at the structure 
of the integrals involved with respect to the variable u. Importantly, the arguments of the Delta functions Si 7^ Si{u) and 
VSi X VSj, given by equation HC12[I . is ~ Qi{u) where Q4,{u) is a positive quartic in u. Puting the expressions for Pi and 
Po into equation (|C9|I imply that the it-component of the integral in dC/dri involves 



I{-1,ri) = 



^Q4(ti)exp(-iiV2) 
(27r)i/2(l_ ^2)1/2 



■ exp 



in - juf 
2(1-72) 



du , 



(C19) 



where Qiiu), of course, also depends on v,w,Uk,Vk,Wk,xi23,Xki with k < I and possibly 7, but not on 7. 

In the trivial limit 7-^0 the coupling between 11 and the field value ry vanishes and the difi'erential length is reduced to 
the PDF of rj 



dC/dr] oc exp 



!L 
2 



Ctot 

(2^ 



■ exp 



■n 



(C20) 



For nonvanishing 7, following NCD, the differentiation of equation (|C19P shows that T{j,t}) obeys the equation 



d_ 

drj 



r?J(7, r?) + 



dl 
drj 



whose solution involve even Hermitte polynomials (retaining only the convergent solution at large rf): 
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n = 



exp 



2 



(C21) 



Note that co is non null since asymptotically, the differential length should converge towards the PDF of 77, within a multi- 
plicative constant. Thanks to the orthogonality condition on the Hermitte polynomial, C2n is given by 



dxH2n{x/V2)exp{~x'^/2)I{x,'^) = / dxH2n{x/V2) exp{ — X^)y/Q4{x) = C2n{v,W,Ui,Ui,Wi,Xi23,Xij) . (C22) 



C2n = lim 

The integration of equation (|C21[I over v,w,Uk,Uk,Wk,xi23,Xki for k < I (while accounting for the rest of the integrant 
corresponding to Po and Pi together with the two Delta functions) yields the functional form of djC/drj, equation (|19|l . where 
C2n is a pure number in the stiff approximation, but may depend on 7 in general. This section falls short of demonstrating 
why only C2 is non null, though the oscillatory behaviour of H2n in equation HC22|) suggests that C2n should decrease with n. 
Making use of the expression of Po and Pi in the stiff regime {xkim ~ 0), equation HC13|) becomes 



dC 
drj 



(1-72)1/2 



exp 



{u - 777)2 
2(1-7^) 



2 



dwexp 



k<l 



3 \ ^ 2 

22^^^ 

k 



d{v, w) 



d(si 



Y[ dxki Y[ dxk 



where v^[a;jfc,Si = 0, Sj = 0] and u?[xjk,Si — 0,Sj = 0] are function of Xki, k < I and Si and Sj (evaluated at zero to 
account for the two Delta functions). The quadratic function Q2{u) corresponds to the square root of Q4{u) when v and w 
are reexpressed in terms of Si (which are in turn evaluated at zero). Now writing formally again (52(1*) — A2U^ + Aiu + ^0, 
if the region where Q2{u) is positive dominates the above integral, the integration over u yields 

[Ao + A2 + Ai-yn + A2 7' H2{v/V2)] exp{-ri^ /2) 

which suggests that C2n ~ for ?i > 1. Given that dC/drj is even in 77, here Ai"f should eventually cancel out when integrated 
over the other variables. 



© 0000 RAS, MNRAS 000, 000-000 



